The data is included within this reproducible report and can be downloaded here (object rbt) and here (object rbt_sci_1).
# create function to invert items
inv7fct <- function(x) (8-as.numeric(x))
# pivot them into long format
rbt_sci_l <- pivot_longer(rbt,
cols = (abs1_tsm_1:abs2_tch_5),
names_to = "variable",
values_to = "value") %>%
mutate(treat = case_when( # create treatment variable
str_detect(variable, "abs1_") ~ toupper(treat1),
str_detect(variable, "abs2_") ~ toupper(treat2)),
variable = str_sub(variable, 6, -1)) # delete title page substring
# pivot them into wide format
rbt_sci_w <- pivot_wider(rbt_sci_l,
names_from = variable,
values_from = value,
id_cols = c(session, treat))
# invert trust items & build scales
data_rbt <- rbt_sci_w %>%
mutate_at(vars(tru_exp_1:tru_ben_4),
list(~inv7fct(.))) %>% # recoding 1-->7, 2-->6, ... %>%
rename(exp_1 = tru_exp_1,
exp_2 = tru_exp_2,
exp_3 = tru_exp_3,
exp_4 = tru_exp_4,
exp_5 = tru_exp_5,
exp_6 = tru_exp_6,
int_1 = tru_int_1,
int_2 = tru_int_2,
int_3 = tru_int_3,
int_4 = tru_int_4,
ben_1 = tru_ben_1,
ben_2 = tru_ben_2,
ben_3 = tru_ben_3,
ben_4 = tru_ben_4) %>%
group_by(session) %>%
mutate(meas_rep = 1:n()) %>%
ungroup() %>%
full_join(.,rbt_sci_1 %>%
select(session, age, sex, position, position_oth) %>%
distinct())
## Joining, by = "session"
# data frames with sum scales
data_scales <- data_rbt%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))
data_scales_lables <- data_rbt%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Treatment = fct_recode(Treatment,
`Colored badges` = "CB",
`Control Condition` = "CC",
`Greyed out badges` = "GB"),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
`Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4),
na.rm = T))
data_scales_wide <- data_scales%>%
select(Experitse, Integrity, Benevolence,
TSM, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
data_scales_lables_wide <- data_scales_lables%>%
select(Experitse, Integrity, Benevolence,
`Topic specific multiplism`, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
# Data of treatment check
data_tch <- data_rbt %>%
select(starts_with("tch"), meas_rep, treat)
data_tch_n <- data_tch%>%
mutate_all(function(x) ifelse(x == -999, NA, x)) %>%
filter(rowSums(is.na(data.frame(tch_1, tch_2, tch_3, tch_4))) != 4)
PID_list <- data_rbt$session %>% unique()
length(PID_list)
## [1] 250
data_rbt %>%
filter(meas_rep == 2) %>%
select(session) %>%
distinct() %>%
nrow()
## [1] 250
sum(is.na(data_rbt%>%
filter(meas_rep == 2)%>%
pull(exp_1)))
## [1] 0
data_rbt %>%
filter(meas_rep == 2)%>%
mutate(position = as.factor(position),
sex = as.factor(sex),
age = as.factor(age)) %>%
select(age, position, sex) %>%
skim(.)
| Name | Piped data |
| Number of rows | 250 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| factor | 3 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| age | 0 | 1 | FALSE | 3 | 16: 193, 35: 37, 50: 20 |
| position | 0 | 1 | FALSE | 6 | -99: 122, 1: 91, 5: 18, 2: 13 |
| sex | 0 | 1 | FALSE | 2 | f: 170, m: 80 |
data_rbt %>%
filter(meas_rep == 2)%>%
pull(sex) %>%
table(.)
## .
## f m
## 170 80
data_rbt %>%
filter(meas_rep == 2)%>%
pull(position) %>%
table(.)
## .
## -999 1 2 3 4 5
## 122 91 13 5 1 18
data_rbt %>%
filter(meas_rep == 2)%>%
pull(age) %>%
table(.)
## .
## 16 35 50
## 193 37 20
skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))
my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))
data_scales%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 500 |
| Number of columns | 38 |
| _______________________ | |
| Column type frequency: | |
| character | 4 |
| factor | 1 |
| numeric | 33 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| session | 0 | 1.00 | 64 | 64 | 0 | 250 | 0 |
| treat | 0 | 1.00 | 2 | 2 | 0 | 3 | 0 |
| sex | 0 | 1.00 | 1 | 1 | 0 | 2 | 0 |
| position_oth | 270 | 0.46 | 2 | 97 | 0 | 103 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Treatment | 0 | 1 | FALSE | 3 | CB: 172, CC: 165, GB: 163 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| tsm_1 | 0 | 1 | 2.03 | 0.83 | 1.00 | 1.00 | 2.00 | 3.00 | 4 | ▆▇▁▅▁ | 0.36 | -0.61 | 2.37 |
| tsm_2 | 0 | 1 | 2.48 | 0.82 | 1.00 | 2.00 | 2.00 | 3.00 | 4 | ▂▇▁▆▂ | 0.24 | -0.51 | 1.85 |
| tsm_3 | 0 | 1 | 2.27 | 0.82 | 1.00 | 2.00 | 2.00 | 3.00 | 4 | ▃▇▁▆▁ | 0.15 | -0.53 | 2.12 |
| tsm_4 | 0 | 1 | 2.67 | 0.80 | 1.00 | 2.00 | 3.00 | 3.00 | 4 | ▁▇▁▇▃ | 0.10 | -0.65 | 2.10 |
| tsc_1 | 0 | 1 | 2.85 | 0.72 | 1.00 | 2.00 | 3.00 | 3.00 | 4 | ▁▅▁▇▂ | -0.15 | -0.32 | 2.56 |
| tsc_2 | 0 | 1 | 1.92 | 0.72 | 1.00 | 1.00 | 2.00 | 2.00 | 4 | ▅▇▁▂▁ | 0.40 | -0.18 | 2.87 |
| tsc_3 | 0 | 1 | 2.90 | 0.70 | 1.00 | 2.75 | 3.00 | 3.00 | 4 | ▁▃▁▇▂ | -0.31 | 0.03 | 2.70 |
| exp_1 | 0 | 1 | 5.49 | 1.29 | 1.00 | 5.00 | 6.00 | 7.00 | 7 | ▁▁▂▃▇ | -0.66 | -0.18 | 3.49 |
| exp_2 | 0 | 1 | 5.56 | 1.17 | 1.00 | 5.00 | 6.00 | 6.00 | 7 | ▁▁▂▃▇ | -0.63 | 0.02 | 3.89 |
| exp_3 | 0 | 1 | 5.66 | 1.18 | 1.00 | 5.00 | 6.00 | 7.00 | 7 | ▁▁▂▃▇ | -0.73 | 0.03 | 3.94 |
| exp_4 | 0 | 1 | 5.47 | 1.37 | 1.00 | 5.00 | 6.00 | 7.00 | 7 | ▁▁▂▃▇ | -0.74 | -0.09 | 3.27 |
| exp_5 | 0 | 1 | 5.13 | 1.38 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▃▅▇ | -0.47 | -0.45 | 3.00 |
| exp_6 | 0 | 1 | 5.40 | 1.28 | 2.00 | 5.00 | 6.00 | 6.00 | 7 | ▂▅▇▇▇ | -0.51 | -0.48 | 2.65 |
| int_1 | 0 | 1 | 5.27 | 1.18 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▃▅▇ | -0.33 | -0.39 | 3.61 |
| int_2 | 0 | 1 | 5.29 | 1.30 | 2.00 | 4.00 | 5.00 | 6.00 | 7 | ▃▆▇▇▆ | -0.41 | -0.58 | 2.53 |
| int_3 | 0 | 1 | 5.10 | 1.20 | 2.00 | 4.00 | 5.00 | 6.00 | 7 | ▂▇▇▇▃ | -0.09 | -0.70 | 2.59 |
| int_4 | 0 | 1 | 5.17 | 1.31 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▅▅▇ | -0.29 | -0.63 | 3.20 |
| ben_1 | 0 | 1 | 5.10 | 1.20 | 2.00 | 4.00 | 5.00 | 6.00 | 7 | ▂▇▇▆▅ | -0.05 | -0.67 | 2.58 |
| ben_2 | 0 | 1 | 5.22 | 1.36 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▃▃▇ | -0.45 | -0.39 | 3.11 |
| ben_3 | 0 | 1 | 5.26 | 1.32 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▃▃▇ | -0.58 | -0.11 | 3.22 |
| ben_4 | 0 | 1 | 5.02 | 1.24 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▆▆▇ | -0.08 | -0.56 | 3.24 |
| tch_1 | 0 | 1 | -111.57 | 318.64 | -999.00 | 1.00 | 2.00 | 4.00 | 4 | ▁▁▁▁▇ | -2.42 | 3.87 | 2.79 |
| tch_2 | 0 | 1 | -215.79 | 413.94 | -999.00 | 1.00 | 1.50 | 4.00 | 4 | ▂▁▁▁▇ | -1.36 | -0.15 | 1.89 |
| tch_3 | 0 | 1 | -209.78 | 409.77 | -999.00 | 1.00 | 1.00 | 4.00 | 4 | ▂▁▁▁▇ | -1.41 | -0.03 | 1.93 |
| tch_4 | 0 | 1 | -133.63 | 343.68 | -999.00 | 1.00 | 2.00 | 4.00 | 4 | ▁▁▁▁▇ | -2.12 | 2.49 | 2.52 |
| tch_5 | 0 | 1 | -211.79 | 411.17 | -999.00 | 1.00 | 1.00 | 4.00 | 4 | ▂▁▁▁▇ | -1.39 | -0.07 | 1.91 |
| meas_rep | 0 | 1 | 1.50 | 0.50 | 1.00 | 1.00 | 1.50 | 2.00 | 2 | ▇▁▁▁▇ | 0.00 | -2.00 | 1.00 |
| age | 0 | 1 | 21.53 | 10.75 | 16.00 | 16.00 | 16.00 | 16.00 | 50 | ▇▁▂▁▁ | 1.67 | 1.36 | 2.65 |
| position | 0 | 1 | -486.61 | 500.74 | -999.00 | -999.00 | 1.00 | 1.00 | 5 | ▇▁▁▁▇ | -0.05 | -2.00 | 1.02 |
| Experitse | 0 | 1 | 5.45 | 1.15 | 1.67 | 4.67 | 5.67 | 6.33 | 7 | ▁▂▅▆▇ | -0.53 | -0.31 | 3.29 |
| Integrity | 0 | 1 | 5.21 | 1.11 | 2.00 | 4.50 | 5.25 | 6.00 | 7 | ▁▅▇▇▆ | -0.24 | -0.58 | 2.90 |
| Benevolence | 0 | 1 | 5.15 | 1.13 | 1.75 | 4.25 | 5.25 | 6.00 | 7 | ▁▃▇▇▇ | -0.23 | -0.48 | 3.01 |
| TSM | 0 | 1 | 2.36 | 0.58 | 1.00 | 2.00 | 2.25 | 2.75 | 4 | ▂▅▇▂▁ | 0.17 | -0.25 | 2.84 |
First we specified two consecutive threedimensional CFA models
# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 2))
cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4"
cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4"
cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_rbt%>%filter(meas_rep == 1))
# define a function which prints the fit
fpf <- function(x){
fm_tmp <- fitmeasures(x)
return(cat(sprintf(
"χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
round(fm_tmp[c("chisq")],3),
fm_tmp[c("df")],
round(fm_tmp[c("cfi")],3),
round(fm_tmp[c("tli")],3),
round(fm_tmp[c("rmsea")],3),
round(fm_tmp[c("srmr")],3),
round(fm_tmp[c("srmr_between")],3),
round(fm_tmp[c("srmr_within")],3))))
}
# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)
χ2 = 463.608, df = 77, CFI = 0.881, TLI = 0.86, RMSEA = 0.142, SRMR = 0.057, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_1)
χ2 = 171.668, df = 74, CFI = 0.97, TLI = 0.963, RMSEA = 0.073, SRMR = 0.029, SRMRbetween = NA, SRMRwithin = NA
cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)
χ2 = 479.509, df = 77, CFI = 0.896, TLI = 0.878, RMSEA = 0.145, SRMR = 0.053, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_2)
χ2 = 166.6, df = 74, CFI = 0.976, TLI = 0.971, RMSEA = 0.071, SRMR = 0.029, SRMRbetween = NA, SRMRwithin = NA
anova(cfa1d_meti_1, cfa_meti_1)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_1 74 8470.3 8579.4 171.67
## cfa1d_meti_1 77 8756.2 8854.8 463.61 291.94 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cfa1d_meti_2, cfa_meti_2)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_2 74 7829.6 7938.8 166.60
## cfa1d_meti_2 77 8136.5 8235.1 479.51 312.91 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In an next step we ran a two-level CFA …
# onedimensional model
mcfa1d_meti_model <- "level: 1
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4
level: 2
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
mcfa_meti_model <- "level: 1
exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_w =~ int_1 + int_2 + int_3 + int_4
ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
level: 2
exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_b =~ int_1 + int_2 + int_3 + int_4
ben_b =~ ben_1 + ben_2 + ben_3 + ben_4
int_2 ~~ int_3
exp_5 ~~ exp_6
int_4 ~~ ben_2
int_3 ~~ ben_1
exp_3 ~~ 0*exp_3
int_3 ~~ 0*int_3
ben_2 ~~ 0*ben_2
ben_3 ~~ 0*ben_3"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_rbt, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_rbt, cluster = "session")
fpf(mcfa1d_meti)
χ2 = 473.759, df = 154, CFI = 0.954, TLI = 0.945, RMSEA = 0.064, SRMR = 0.304, SRMRbetween = 0.232, SRMRwithin = 0.072
fpf(mcfa_meti)
χ2 = 238.803, df = 148, CFI = 0.987, TLI = 0.984, RMSEA = 0.035, SRMR = 0.125, SRMRbetween = 0.085, SRMRwithin = 0.04
anova(mcfa1d_meti, mcfa_meti)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## mcfa_meti 148 16218 16538 238.80
## mcfa1d_meti 154 16441 16736 473.76 234.96 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# First Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("exp")))$est
## [1] 0.9420418
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("int")))$est
## [1] 0.9084106
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("ben")))$est
## [1] 0.9006582
# Second Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("exp")))$est
## [1] 0.9623451
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("int")))$est
## [1] 0.9179565
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("ben")))$est
## [1] 0.9110896
First we specified two consecutive onedimensional CFA models
cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
tsm_2 ~~ tsm_4"
cfa_tsm_1 <-
cfa(cfa_tsm_model,
data = data_rbt %>%
filter(meas_rep == 1),
missing = "fiml")
fpf(cfa_tsm_1)
χ2 = 5.245, df = 4, CFI = 0.993, TLI = 0.99, RMSEA = 0.035, SRMR = 0.034, SRMRbetween = NA, SRMRwithin = NA
cfa_tsm_2 <-
cfa(cfa_tsm_model,
data = data_rbt %>%
filter(meas_rep == 2))
fpf(cfa_tsm_2)
χ2 = 7.677, df = 4, CFI = 0.974, TLI = 0.962, RMSEA = 0.061, SRMR = 0.058, SRMRbetween = NA, SRMRwithin = NA
In an next step, we ran a two-level CFA …
mcfa_tsm_model <- "level: 1
tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ tsm_4
level: 2
tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ tsm_4
"
mcfa_tsm <- cfa(mcfa_tsm_model, data = data_rbt, cluster = "session")
fpf(mcfa_tsm)
χ2 = 1.023, df = 2, CFI = 1, TLI = 1.02, RMSEA = 0, SRMR = 0.027, SRMRbetween = 0.019, SRMRwithin = 0.008
# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==1) %>% select(starts_with("tsm")))$est
## [1] 0.6911811
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==2) %>% select(starts_with("tsm")))$est
## [1] 0.6447776
We specified two consecutive onedimensional CFA models
cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5
tch_1 ~~ tch_4"
cfa_tch_model2 <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5
tch_1 ~~ tch_4
tch_3 ~~ tch_5
tch_1 ~~ tch_2"
cfa_tch_1 <- cfa(cfa_tch_model,
data = data_tch_n%>%
filter(meas_rep == 1))
fpf(cfa_tch_1)
χ2 = 9.287, df = 4, CFI = 0.996, TLI = 0.989, RMSEA = 0.089, SRMR = 0.008, SRMRbetween = NA, SRMRwithin = NA
cfa_tch_2 <- cfa(cfa_tch_model2, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)
χ2 = 0.269, df = 2, CFI = 1, TLI = 1.006, RMSEA = 0, SRMR = 0.001, SRMRbetween = NA, SRMRwithin = NA
modificationindices(cfa_tch_2, sort. = T)
lhs op rhs mi epc sepc.lv sepc.all sepc.nox
20 tch_3 ~~ tch_4 0.262 0.009 0.009 0.045 0.045 21 tch_4 ~~ tch_5 0.262 -0.009 -0.009 -0.050 -0.050 19 tch_2 ~~ tch_5 0.075 0.005 0.005 0.043 0.043 17 tch_2 ~~ tch_3 0.075 -0.005 -0.005 -0.038 -0.038 15 tch_1 ~~ tch_3 0.072 -0.005 -0.005 -0.020 -0.020 16 tch_1 ~~ tch_5 0.072 0.005 0.005 0.022 0.022
# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 1) %>% select(starts_with("tch")))$est
## [1] 0.8746516
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est
## [1] 0.9149272
tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
column_to_rownames(var = "rownames")%>%
knitr::kable(., digits = 3)
| 1d CFA METI 1 | 1d CFA METI 2 | 3d CFA METI 1 | 3d CFA METI 2 | 1d MCFA METI | 3d MCFA METI | 1d CFA TSM 1 | 1d CFA TSM 2 | 1d MCFA TSM | 1d CFA TCH 1 | 1d CFA TCH 2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| χ2 | 463.608 | 479.509 | 171.668 | 166.600 | 473.759 | 238.803 | 5.245 | 7.677 | 1.023 | 9.287 | 0.269 |
| df | 77.000 | 77.000 | 74.000 | 74.000 | 154.000 | 148.000 | 4.000 | 4.000 | 2.000 | 4.000 | 2.000 |
| CFI | 0.881 | 0.896 | 0.970 | 0.976 | 0.954 | 0.987 | 0.993 | 0.974 | 1.000 | 0.996 | 1.000 |
| TLI | 0.860 | 0.878 | 0.963 | 0.971 | 0.945 | 0.984 | 0.990 | 0.962 | 1.020 | 0.989 | 1.006 |
| RMSEA | 0.142 | 0.145 | 0.073 | 0.071 | 0.064 | 0.035 | 0.035 | 0.061 | 0.000 | 0.089 | 0.000 |
| SRMR | 0.057 | 0.053 | 0.029 | 0.029 | 0.304 | 0.125 | 0.034 | 0.058 | 0.027 | 0.008 | 0.001 |
| SRMRbetween | NA | NA | NA | NA | 0.232 | 0.085 | NA | NA | 0.019 | NA | NA |
| SRMRwithin | NA | NA | NA | NA | 0.072 | 0.040 | NA | NA | 0.008 | NA | NA |
| BIC | 8854.815 | 8235.105 | 8579.439 | 7938.761 | 16735.574 | 16537.906 | 2369.599 | 2237.418 | 4559.431 | 1667.935 | 1685.730 |
| AIC | 8756.214 | 8136.504 | 8470.274 | 7829.596 | 16440.552 | 16217.596 | 2334.384 | 2216.290 | 4466.709 | 1633.572 | 1644.514 |
res_tch_data <- data_tch%>%
gather(variable, value, starts_with("tch_"))%>%
group_by(treat, variable, value)%>%
mutate(value = as.character(value)) %>%
summarize(freq = n())%>%
ungroup()%>%
mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ treat),
value = case_when(value =="-999" ~ "don't know",
T ~ value),
variable = case_when(variable == "tch_1" ~ "item 1",
variable == "tch_2" ~ "item 2",
variable == "tch_3" ~ "item 3",
variable == "tch_4" ~ "item 4",
variable == "tch_5" ~ "item 5"),
Frequency = freq)
## `summarise()` regrouping output by 'treat', 'variable' (override with `.groups` argument)
res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) +
geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
scale_size_continuous(range = c(3,15)) +
scale_color_gradient(low = "grey95", high = "grey65") +
guides(color=guide_legend(), size = guide_legend()) +
facet_wrap(~treat) +
theme_ipsum_ps() +
ggtitle("Results of the treatment check", "Frequency per item and experimental condition") +
ylab("") +
xlab("")
res_tch_plot
res_tch_plot_pub <- ggplot(res_tch_data %>%
mutate(treat = case_when(treat == "Greyed out badges" ~ "Grey badges",
TRUE ~ treat)), aes(variable, value)) +
geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
scale_size_continuous(range = c(3,15)) +
scale_color_gradient(low = "grey95", high = "grey65") +
guides(color=guide_legend(), size = guide_legend()) +
facet_wrap(~treat) +
theme_ipsum_ps() +
ylab("") +
xlab("")
#ggsave("8_reproducible_documentation_of_analyses/Fig/Fig6.jpg", res_tch_plot_pub, width = 130, height = 50, units = "mm", dpi = 300, scale = 2.4)
res_tch_data_A <- data_tch_n%>%
filter(treat != "CC")%>%
filter(tch_1 != -999)
effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)
##
## Vargha and Delaney A
##
## A estimate: 0.9378639 (large)
data_scales_lables%>%
gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
ggplot(., aes(x = Treatment, y = Value)) +
geom_violin(adjust = 1.5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
facet_wrap(~ Variable, nrow = 1) +
labs(title = "Graphical overview (Hyp 1)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
ylim(1,7) +
hrbrthemes::theme_ipsum_ps()
# Export data for combined plot
#data_scales_lables %>%
# mutate(Study = "Study 2") %>%
# write_csv(., "8_reproducible_documentation_of_analyses/Fig/data_scales_lables_study_2.csv")
A_GB_CC <- data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_CC_CB <- data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_GB_CB <- data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
| GB | CC | |
|---|---|---|
| CC | A = 0.65, d = 0.55 | |
| CB | A = 0.71, d = 0.77 | A = 0.57, d = 0.25 |
plot_hyp2_1 <- ggplot(data_scales_lables,
aes(x=`Topic specific multiplism`, y = Integrity)) +
geom_jitter() +
facet_wrap(~ Treatment, nrow = 1) +
labs(title = "Graphical overview (Hyp 2/3)",
subtitle = "Jitter plot per treatment") +
hrbrthemes::theme_ipsum_ps()
plot_hyp2_1 + stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
plot_hyp2_1 + stat_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'
Spearman and Kendall correlations:
r_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
| GB | CC | CB | |
|---|---|---|---|
| \(r(integrity, multiplism)\) | -0.34 | -0.32 | -0.35 |
| \(\tau(integrity, multiplism)\) | -0.23 | -0.23 | -0.24 |
data_scales_lables%>%
mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ "treat"))%>%
ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
geom_violin(adjust = 1.5, alpha = .5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
labs(title = "Graphical overview (Hyp 4)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
xlab("") +
ylim(1,4) +
hrbrthemes::theme_ipsum_ps()
fig4 <- data_scales_lables%>%
mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ "treat"))%>%
ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
geom_violin(adjust = 1.5, alpha = .5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
xlab("") +
ylim(1,4) +
hrbrthemes::theme_ipsum_ps()
A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)
A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)
A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)
| GB | CC | |
|---|---|---|
| CC | A = 0.43, d = -0.27 | |
| CB | A = 0.43, d = -0.24 | A = 0.51, d = 0.02 |
All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:
GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open PracticesCC means, that the participants of the study could not learn if or if not the authors of the study used Open PracticesCBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open PracticesFinally, the variable session identified the study participants.
If we look descriptively at these variables:
data_scales_wide%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 250 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| session | 0 | 1 | 64 | 64 | 0 | 250 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Benevolence_CB | 78 | 0.69 | 5.51 | 1.05 | 2.75 | 4.75 | 5.50 | 6.25 | 7 | ▁▃▇▇▇ | -0.34 | -0.63 | 2.63 |
| Benevolence_CC | 85 | 0.66 | 5.27 | 1.02 | 2.50 | 4.50 | 5.25 | 6.00 | 7 | ▁▅▆▇▅ | -0.10 | -0.68 | 2.73 |
| Benevolence_GB | 87 | 0.65 | 4.64 | 1.15 | 1.75 | 3.75 | 4.75 | 5.50 | 7 | ▂▆▇▇▅ | -0.06 | -0.43 | 2.52 |
| Experitse_CB | 78 | 0.69 | 5.72 | 1.05 | 2.33 | 5.12 | 5.83 | 6.67 | 7 | ▁▂▃▇▇ | -0.65 | -0.22 | 3.23 |
| Experitse_CC | 85 | 0.66 | 5.71 | 1.04 | 2.50 | 5.17 | 5.83 | 6.50 | 7 | ▁▂▃▇▇ | -0.80 | 0.40 | 3.09 |
| Experitse_GB | 87 | 0.65 | 4.91 | 1.18 | 1.67 | 4.00 | 5.00 | 5.83 | 7 | ▁▃▇▇▆ | -0.11 | -0.49 | 2.76 |
| Integrity_CB | 78 | 0.69 | 5.57 | 1.01 | 3.00 | 4.75 | 5.75 | 6.50 | 7 | ▁▅▅▇▇ | -0.27 | -0.93 | 2.55 |
| Integrity_CC | 85 | 0.66 | 5.32 | 0.99 | 2.75 | 4.50 | 5.25 | 6.00 | 7 | ▁▅▇▇▆ | -0.12 | -0.84 | 2.61 |
| Integrity_GB | 87 | 0.65 | 4.71 | 1.15 | 2.00 | 4.00 | 4.75 | 5.50 | 7 | ▂▆▇▆▃ | -0.03 | -0.54 | 2.36 |
| TSM_CB | 78 | 0.69 | 2.31 | 0.59 | 1.00 | 2.00 | 2.25 | 2.75 | 4 | ▂▅▇▂▁ | -0.01 | -0.42 | 2.86 |
| TSM_CC | 85 | 0.66 | 2.32 | 0.59 | 1.00 | 2.00 | 2.25 | 2.75 | 4 | ▂▆▇▂▁ | 0.37 | -0.17 | 2.87 |
| TSM_GB | 87 | 0.65 | 2.46 | 0.54 | 1.25 | 2.00 | 2.50 | 2.75 | 4 | ▃▇▇▅▁ | 0.26 | -0.30 | 2.83 |
library(naniar)
##
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
##
## n_complete
visdat::vis_miss(data_scales_wide%>%select(-session)) +
ggtitle("Missingness per Variable") +
theme(plot.margin = margin(1, 2, 1, 1, "cm"))
Integrity_GBlibrary(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CCggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CBggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_GBggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CCggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CBggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
round(.,4)
knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")
| Integrity_CB | Integrity_CC | Integrity_GB | TSM_CB | TSM_CC | TSM_GB | |
|---|---|---|---|---|---|---|
| Integrity_CB | NaN | -0.0349 | -0.0827 | NaN | -0.2 | -0.0511 |
| Integrity_CC | 0.0504 | NaN | 0.0827 | -0.0026 | NaN | 0.0511 |
| Integrity_GB | -0.0504 | 0.0349 | NaN | 0.0026 | 0.2 | NaN |
| TSM_CB | NaN | -0.0349 | -0.0827 | NaN | -0.2 | -0.0511 |
| TSM_CC | 0.0504 | NaN | 0.0827 | -0.0026 | NaN | 0.0511 |
| TSM_GB | -0.0504 | 0.0349 | NaN | 0.0026 | 0.2 | NaN |
boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")
M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
m = M,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
out_first10 <- mice(data = data_scales_wide%>%select(-session),
m = 10,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
densityplot(out_first10)
plot(out_first10)
# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3)
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3)
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance
# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))*
sum(diag(covbetween_hyp1 %*%
MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)
# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp1 <- bain(estimates_hyp1,
"Integrity_GB<Integrity_CC<Integrity_CB;
Integrity_GB=Integrity_CC=Integrity_CB;
Integrity_GB<Integrity_CC=Integrity_CB",
n = neff_hyp1, Sigma=covariance_hyp1,
group_parameters=3, joint_parameters = 0)
print(results_hyp1)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.994 0.183 5.438 803.908 0.882 0.759
## H2 0.000 0.203 0.000 0.000 0.000 0.000
## H3 0.166 0.229 0.725 0.725 0.118 0.101
## Hu 0.140
##
## Hypotheses:
## H1: Integrity_GB<Integrity_CC<Integrity_CB
## H2: Integrity_GB=Integrity_CC=Integrity_CB
## H3: Integrity_GB<Integrity_CC=Integrity_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
summary(results_hyp1)
## Parameter n Estimate lb ub
## 1 Integrity_GB 167.6928 4.728577 4.556397 4.900757
## 2 Integrity_CC 167.6928 5.319534 5.174336 5.464732
## 3 Integrity_CB 167.6928 5.561620 5.413517 5.709722
results_hyp1$BFmatrix
## H1 H2 H3
## H1 1.000000e+00 155466504142 7.497718e+00
## H2 6.432254e-12 1 4.822723e-11
## H3 1.333739e-01 20735175670 1.000000e+00
path_mod <- "Integrity_GB ~ TSM_GB
Integrity_CC ~ TSM_CC
Integrity_CB ~ TSM_CB"
# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
path_fit <- sem(path_mod,
data = mice::complete(out, i),
std.lv = T) # estimate the path coefficients
best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
filter(op == "~")%>%
pull(std.all)
covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
1/M * lavInspect(path_fit,
"vcov.std.all")[c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB"),
c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB")]
}
# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)
## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB
## -0.35 -0.37 -0.41
\usetikzlibrary{arrows,positioning}
\usetikzlibrary{calc}
\begin{tikzpicture}[auto,>=latex,align=center,
latent/.style={circle,draw, thick,inner sep=0pt,minimum size=10mm},
manifest/.style={rectangle,draw, thick,inner sep=2pt,minimum size=15mm},
manifestcov/.style={rectangle,draw, thick,inner sep=1pt,minimum size=8mm},
mean/.style={regular polygon,regular polygon sides=3,draw,thick,inner sep=0pt,minimum size=10mm},
paths/.style={->, thick, >=latex, font = \footnotesize\sffamily, fill = white},
variance/.style={<->, thin, >=latex, bend left=270, looseness=2.5, font = \footnotesize},
covariance/.style={<->, thin, >=latex, bend left=290, looseness=.5, font = \footnotesize},
dotsdots/.style={circle, draw, fill = black, minimum size = 1mm, maximum size = 1mm}
]
%% AVs
\node [manifest] (AV1) at (0, -2) {$int_{gb}$};
\node [manifest] (AV2) [right = of AV1] {$int_{cc}$};
\node [manifest] (AV3) [right = of AV2] {$int_{cb}$};
%% UVs
\node [manifest] (UV1) at (0, 2) {$tsm_{gb}$};
\node [manifest] (UV2) [right =of UV1] {$tsm_{cc}$};
\node [manifest] (UV3) [right =of UV2] {$tsm_{cb}$};
%% Paths
\draw [paths] (UV1) to node [midway]{-.35} (AV1);
\draw [paths] (UV2) to node [midway]{-.37} (AV2);
\draw [paths] (UV3) to node [midway]{-.41} (AV3);
%\draw [covariance] (UV2) to node {} (UV1);
\end
{tikzpicture}
covbetween_hyp2 <- cov(best_hyp2) # between covariance matrix
covariance_hyp2 <- covwithin_hyp2 + (1+1/M)*covbetween_hyp2 # total variance
# Determine the effective and real sample sizes ######
samp_hyp2 <- nrow(data_scales_wide) # real sample size
nucom_hyp2 <- samp_hyp2-length(estimates_hyp2)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp2 <- (1+1/M)*(1/length(estimates_hyp2))*
sum(diag(covbetween_hyp2 %*%
MASS::ginv(covariance_hyp2))) # ... (43)
nuold_hyp2 <- (M-1)/(lam_hyp2^2) # ... (44)
nuobs_hyp2 <- (nucom_hyp2+1)/(nucom_hyp2+3)*nucom_hyp2*(1-lam_hyp2) # ... (46)
nu_hyp2 <- nuold_hyp2*nuobs_hyp2/(nuold_hyp2+nuobs_hyp2) # ... (47)
fracmis_hyp2 <- (nu_hyp2+1)/(nu_hyp2+3)*lam_hyp2 + 2/(nu_hyp2+3) # ... (48)
neff_hyp2 <- samp_hyp2-samp_hyp2*fracmis_hyp2 # = 114 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp2 <- list(covariance_hyp2)
set.seed(6346)
results_hyp2 <- bain(estimates_hyp2,
"Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = F)
print(results_hyp2)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 1.000 0.145 6.890 26134059.789 1.000 0.873
## H2 0.000 0.741 0.000 0.000 0.000 0.000
## Hu 0.127
##
## Hypotheses:
## H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
## H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp2$BFmatrix
## H1 H2
## H1 1.000000e+00 2.558656e+16
## H2 3.908302e-17 1.000000e+00
set.seed(6346)
results_hyp3 <- bain(estimates_hyp2,
"Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
(Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = T)
print(results_hyp3)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.069 0.164 0.420 0.377 0.008 0.008
## H2 0.145 0.322 0.451 0.358 0.008 0.008
## H3 18.337 0.341 53.716 53.716 0.984 0.966
## Hu 0.018
##
## Hypotheses:
## H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
## H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
## H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp3$BFmatrix
## H1 H2 H3
## H1 1.00000 0.9313845 0.007820809
## H2 1.07367 1.0000000 0.008396971
## H3 127.86401 119.0905634 1.000000000
# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance
# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))*
sum(diag(covbetween_hyp4 %*%
MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)
# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp4 <- bain(estimates_hyp4,
"TSM_GB>TSM_CC>TSM_CB;
TSM_GB=TSM_CC=TSM_CB;
(TSM_GB,TSM_CC)>TSM_CB",
n = neff_hyp4, Sigma=covariance_hyp4,
group_parameters=3, joint_parameters = 0)
print(results_hyp4)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.742 0.161 4.621 15.055 0.648 0.568
## H2 0.156 0.755 0.206 0.206 0.029 0.025
## H3 0.751 0.325 2.308 6.252 0.323 0.284
## Hu 0.123
##
## Hypotheses:
## H1: TSM_GB>TSM_CC>TSM_CB
## H2: TSM_GB=TSM_CC=TSM_CB
## H3: (TSM_GB,TSM_CC)>TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp4$BFmatrix
## H1 H2 H3
## H1 1.00000000 22.41071 2.00252396
## H2 0.04462152 1.00000 0.08935566
## H3 0.49936981 11.19123 1.00000000